sc10x analysis
load dependancies
Baf53cre Ahr-cko vs ctl,
CR infection,
Intestinal RORgt+ Immune Cells
loading 72,000 for all
cellranger calling 34,620 cells, mean reads 11,163
t1 advanced
load obj
GEX.seur <- readRDS("./20230927_10x_SZJ.new20270809.preAnno.rds")
GEX.seur## An object of class Seurat
## 16621 features across 24449 samples within 1 assay
## Active assay: RNA (16621 features, 1224 variable features)
## 3 dimensional reductions calculated: pca, tsne, umap#
color.FB <- ggsci::pal_d3("category20c")(20)[c(2,7,12,17,
1,6,11,16
)]
color.cnt <- color.FB[c(5,1)]DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno", label.size = 2.75) +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
filtering
# remove Mix
# remove DF0.1-Doublets except cycling ILC3/Th
GEX.seur <- subset(GEX.seur, subset = preAnno != "Mix")
GEX.seur <- subset(GEX.seur, subset = DoubletFinder0.1=="Singlet" | preAnno %in% c("Th.CC","ILC3.CC"))
GEX.seur## An object of class Seurat
## 16621 features across 22088 samples within 1 assay
## Active assay: RNA (16621 features, 1224 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapDimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno", label.size = 2.75) +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
VlnPlot(subset(GEX.seur, subset = cnt %in% c("CTL")),
features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "preAnno")
VlnPlot(subset(GEX.seur, subset = cnt %in% c("CKO")),
features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "preAnno")
markers.new <- c("Cd3d","Cd3e","Cd4","Cd8a",
"Klf2","Ccr7","Lef1","S1pr1",
"Tubb2a","Rflnb","Igfbp4","Sell",
"Itm2a","Itga6","Smc4",#"Sox4","Pdcd1",
"Ptprc",#"Fam241a","Slamf6",
"Izumo1r","Tbc1d4",
"Ctla4","Il10",
"Hopx","Foxp3","Ccr2","Maf",
"Capg","Slamf1","Icos",
"Il17a","Nebl","S100a11","S100a10",
"Top2a","Mki67","Pcna","Mcm6",
"Fasl","Gzma","Gzmk","Ccr5",
"Il12rb2",
"Jun","Fos","Ier2","Egr1",
"Ifng", "Tnf",
"Tbx21",
"Klrb1c","Xcl1","Nkg7","Plac8","Cd160",
"Ccl5",
"Klrd1","Cd7","Itgae",
"Tcrg-C1","Trdv4","Cd163l1","Blk","Mmp25",
"Calca","Gata3","Il17rb","Il13",
"Csf2","Dgat2",
"Il4","Hilpda","Ar","Il1rl1",
"Il5",#"Cd24a",
"Ccl1","Areg","Ahr",
"Cxcl10",
"Gbp6","Stat1","Isg15","Ifit1",
"Ncr1","Gzmc","Gzmb","Irf8",
"Cxcr6","Irf7","Slc16a6","Klrk1",
"Tmem176b","Fcer1g","Car2","Ckb",
"Cd69","Jaml","Cxcl9",
"Il22","Gem","Pxdc1","Cdc14a",
"Sdc4","Vegfa","Bst2","Nmrk1",
"Hmgn3","Cd81","Cx3cl1",
"Cd74","H2-Aa","H2-Eb1","H2-Ab1",
"Ccr6","Batf3","Rorc","Maff",
"Icam1","Il18r1",
"Il17f"
)
DotPlot(GEX.seur, features = markers.new, group.by = "preAnno",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) +
scale_y_discrete(limits=rev)
re-clustering
GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)
# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2
head(VariableFeatures(GEX.seur), 200)## [1] "Hist1h1b" "Il22" "Ccl4" "Il17a" "Gzma" "Ccl1"
## [7] "Ccl5" "Cxcl2" "Cxcl10" "Areg" "Xcl1" "Hist1h2ae"
## [13] "Calca" "Stmn1" "Cxcl1" "Hspa1a" "S100g" "Ube2c"
## [19] "Il13" "Il5" "Ccl20" "Gzmc" "Il4" "Cd74"
## [25] "Pclaf" "S100a6" "Gzmb" "Cxcl3" "Cxcl9" "Tnfrsf4"
## [31] "Il17f" "Ifng" "Lgals1" "Ccl3" "Ifi27l2a" "Hs3st1"
## [37] "Il10" "Trdc" "Egr1" "Penk" "H2-Aa" "Top2a"
## [43] "Mki67" "Stmn2" "Csf2" "Il2" "Pcp4" "Hist1h2ap"
## [49] "Tph1" "Birc5" "Plac8" "Defa24" "H2-Eb1" "Hmgb2"
## [55] "Hist1h3c" "Cenpf" "Serpine1" "Sostdc1" "Hspb1" "Ctla2a"
## [61] "Jun" "Egr2" "Hist1h2ab" "Serpinb1a" "Lgals3" "Ifitm2"
## [67] "Rrm2" "Try5" "S100a4" "Ifitm1" "Bambi" "Cd36"
## [73] "Tmsb4x" "Rrad" "Ccr7" "Ctla4" "Actb" "Crip1"
## [79] "Tubb5" "Hspa1b" "Ifit1" "Ly6a" "Hist1h1e" "Gzmf"
## [85] "Bcl2a1b" "Cdca8" "Fabp5" "Gzme" "Nusap1" "H2-Ab1"
## [91] "Tubb2a" "Akap12" "Klf2" "Vim" "G0s2" "Mt1"
## [97] "Dppa2" "Nkg7" "Rgcc" "Defa21" "Tnfsf8" "Igkc"
## [103] "Fos" "Apoe" "Spp1" "C2cd4b" "Cd7" "Ptgs2"
## [109] "Tyrobp" "Egr3" "Ifitm3" "Lyz1" "Hbegf" "Igha"
## [115] "Ccnb2" "Dnaja1" "Cldn5" "Hes1" "Cd83" "Gadd45g"
## [121] "Ms4a4b" "Cks1b" "Gata3" "Actg2" "Rasl11a" "Gm10827"
## [127] "Defa22" "Gadd45b" "Isg15" "Gfod2" "Tyms" "Clspn"
## [133] "Jchain" "Trac" "Gzmd" "Trbc1" "Klra5" "Iigp1"
## [139] "Hopx" "Tuba1b" "Ptpn13" "Arg1" "Hmmr" "Hba-a1"
## [145] "Il21" "Defa30" "Ermn" "Fbxo5" "Ccna2" "Dennd5b"
## [151] "Hilpda" "Klra1" "Muc3" "Atf3" "Spc24" "Odc1"
## [157] "Klrd1" "AY761184" "Il1r2" "Slpi" "Dnajb1" "Il1rl1"
## [163] "Cdca3" "Acod1" "Cd3g" "Fst" "Tpx2" "Klra7"
## [169] "H2afz" "Dut" "Bcl2a1a" "Kif11" "Hist1h1d" "Cd5"
## [175] "Il6" "Cd24a" "Zfp36" "Gm4841" "Serpinb9" "Igfbp4"
## [181] "Fabp4" "Cd163l1" "Trdv4" "Esco2" "Socs2" "Cenpe"
## [187] "Egr4" "Chad" "Eprs" "Gm14851" "Serpinb6b" "Smc2"
## [193] "Id2" "Izumo1r" "Arl5c" "Gem" "Tent5a" "Ccr2"
## [199] "Akr1c18" "Cd3d"GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))## Centering and scaling data matrix# exclude MT genes and more
DIG <- grep("^mt-|^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^AC|^AA|^AW|^BC|^Gm|^Hist|Rik$|-ps|^Ifi|^Isg|^Egr",
# try adding some few to reduce the mixed Ifit+
rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
c(#MT_gene,
DIG,
CC_gene) )
GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))## PC_ 1
## Positive: Fcer1g, Tmem176a, Gem, Tmem176b, Itm2b, Pxdc1, Sdc4, Il22, Car2, Nrgn
## Gpr183, Cdkn1a, Eprs, Zfp36l1, Ckb, Prr29, Ccdc184, Klrb1b, Hk2, Odc1
## Klf4, Cwc25, Tnfsf11, Rgs1, Atf3, Fam110a, Arrdc4, Gpx1, Dad1, Ostf1
## Negative: Cd3d, Cd3g, Cd2, Gramd3, Ms4a4b, Fyb, Gimap6, Cd28, Cd3e, Ms4a6b
## Il21r, Dusp10, Lef1, Klf2, Lat, Ccr7, Pmepa1, Slfn2, Satb1, S1pr1
## Stk4, Gimap4, Gimap1, Cd5, Cd27, Coro1a, Cd247, Sh2d2a, Got1, Gimap7
## PC_ 2
## Positive: Ccr7, Vegfa, Cdc14a, Hmgn3, Cd81, Bst2, S1pr1, Sdc4, Ramp1, Cx3cl1
## Klf2, Cd74, Batf3, Igfbp4, Lef1, Cryaa, Rflnb, Bmp2, Irs2, Rgcc
## Myof, Nrp1, Ccr6, Cited2, Ttn, H2-Aa, Lmna, Cd83, Il22, Tgm2
## Negative: Tmsb4x, Bcl2a1d, Cxcr6, Serpinb9, Pfn1, Lgals1, Capg, Serpinb6b, Ucp2, Pik3r1
## Icos, Ctsw, Cd52, Serpinb1a, Bcl2a1b, Actn2, Srgn, Ltb4r1, Slc16a6, Ikzf2
## Ptprcap, Maf, Clic1, Slc7a6os, S100a6, Rgs1, Nkg7, S100a11, Rbpms2, Arl5c
## PC_ 3
## Positive: Hs3st1, Calca, Ptpn13, Il5, Gata3, Ccl1, Klrg1, Il17rb, Il13, Il1rl1
## Ipmk, Csf2, Areg, Il4, Slc7a8, Otulinl, Tspan13, Deptor, Lmo4, Klf5
## Tnfrsf18, Hes1, Spcs3, Rab27b, Tent5a, Mras, Rab4a, Ptgir, Ar, Dgat2
## Negative: Fcer1g, Ckb, Ms4a4b, Lef1, Car2, Ccr7, Serpinb1a, Gpx1, Ms4a6b, Gimap4
## S1pr1, Klf2, Tubb2a, Klrk1, Gramd3, Capg, Cd3d, Serpinb6b, Serpinb9, Rflnb
## Prr29, Stmn1, Klrb1b, Il22, Upp1, Calm1, Rbpms2, Cd2, Rgcc, Slc7a6os
## PC_ 4
## Positive: Pclaf, Spc24, Actb, Stmn1, Tnfrsf4, Ccna2, Asf1b, Tk1, Racgap1, Kif4
## Vim, S100a4, S100a6, Ccr6, Ctla4, Knl1, Smc2, Actg1, Lgals1, Bst2
## Mcm3, Cenpm, Lig1, Hmgn3, Cd74, Dut, Ramp1, Il17a, Ly6a, Tubb5
## Negative: Serpinb9, Serpinb6b, Pik3r1, Ctsw, Lef1, Upp1, Klrk1, Satb1, Slc7a6os, Rbpms2
## Tubb2a, Ncr1, Rflnb, Gimap6, Gimap4, Irf7, Dusp10, Ikzf2, Pcp4, Slc16a10
## Slc16a6, Serpinb1a, Klf2, Sell, Lmo4, Ccr7, Gjb2, Ppfia1, Cd200r4, Trf
## PC_ 5
## Positive: Bcl2a1b, Cd3g, Ccr2, Ctla4, Cd3e, Hopx, Ccl5, S100a6, Tnfrsf4, S100a10
## Ccr5, Got1, Fth1, Cd6, Actb, Il10, Odc1, Cd52, Actn2, Cd163l1
## S100a11, Cd3d, Foxp3, Sh2d2a, Plac8, Ltb, Blk, Ubald2, Nkg7, Cd40lg
## Negative: Pclaf, Stmn1, Spc24, Ccna2, Asf1b, Kif4, Knl1, Cenpm, Tk1, Smc2
## Shcbp1, Spc25, Tubb5, Lef1, Racgap1, Ccr7, Nrm, Kif15, S1pr1, Rflnb
## Cdkn3, Tubb2a, Pbk, Tuba1b, Cenph, Aspm, Mcm3, Hmgn2, Esco2, Hs3st1length(setdiff(VariableFeatures(object = GEX.seur),
c(#MT_gene,
DIG,CC_gene)))## [1] 1197setdiff(VariableFeatures(object = GEX.seur),
c(#MT_gene,
DIG,CC_gene))[1:300]## [1] "Il22" "Ccl4" "Il17a" "Gzma" "Ccl1" "Ccl5"
## [7] "Cxcl2" "Cxcl10" "Areg" "Xcl1" "Calca" "Stmn1"
## [13] "Cxcl1" "S100g" "Il13" "Il5" "Ccl20" "Gzmc"
## [19] "Il4" "Cd74" "Pclaf" "S100a6" "Gzmb" "Cxcl3"
## [25] "Cxcl9" "Tnfrsf4" "Il17f" "Ifng" "Lgals1" "Ccl3"
## [31] "Hs3st1" "Il10" "Penk" "H2-Aa" "Stmn2" "Csf2"
## [37] "Il2" "Pcp4" "Tph1" "Plac8" "Defa24" "H2-Eb1"
## [43] "Serpine1" "Sostdc1" "Ctla2a" "Serpinb1a" "Lgals3" "Try5"
## [49] "S100a4" "Bambi" "Cd36" "Tmsb4x" "Rrad" "Ccr7"
## [55] "Ctla4" "Actb" "Crip1" "Tubb5" "Ly6a" "Gzmf"
## [61] "Bcl2a1b" "Fabp5" "Gzme" "H2-Ab1" "Tubb2a" "Akap12"
## [67] "Klf2" "Vim" "G0s2" "Mt1" "Dppa2" "Nkg7"
## [73] "Rgcc" "Defa21" "Tnfsf8" "Apoe" "Spp1" "C2cd4b"
## [79] "Cd7" "Ptgs2" "Tyrobp" "Lyz1" "Hbegf" "Cldn5"
## [85] "Hes1" "Cd83" "Gadd45g" "Ms4a4b" "Gata3" "Actg2"
## [91] "Rasl11a" "Defa22" "Gadd45b" "Gfod2" "Gzmd" "Klra5"
## [97] "Iigp1" "Hopx" "Tuba1b" "Ptpn13" "Arg1" "Il21"
## [103] "Defa30" "Ermn" "Fbxo5" "Ccna2" "Dennd5b" "Hilpda"
## [109] "Klra1" "Muc3" "Atf3" "Spc24" "Odc1" "Klrd1"
## [115] "Il1r2" "Slpi" "Il1rl1" "Acod1" "Cd3g" "Fst"
## [121] "Klra7" "H2afz" "Dut" "Bcl2a1a" "Cd5" "Il6"
## [127] "Cd24a" "Zfp36" "Serpinb9" "Igfbp4" "Fabp4" "Cd163l1"
## [133] "Esco2" "Socs2" "Chad" "Eprs" "Serpinb6b" "Smc2"
## [139] "Id2" "Izumo1r" "Arl5c" "Gem" "Tent5a" "Ccr2"
## [145] "Akr1c18" "Cd3d" "Lef1" "Sptssb" "Gramd3" "Blnk"
## [151] "Nrgn" "Lyz2" "Tnfsf11" "Lgals7" "Klra6" "Lmo4"
## [157] "Ier3" "Nlrp3" "Adra2a" "Bcl2a1d" "Cd2" "Upp1"
## [163] "Vegfa" "Cd52" "Adm" "Prss2" "Tpm2" "Cd28"
## [169] "Tnfrsf9" "Dusp14" "Basp1" "Rsad2" "Serpine2" "Rflnb"
## [175] "Actn2" "Gcg" "Msx1" "Mcm3" "Ptprz1" "Cst8"
## [181] "Il1rn" "Evx1os" "Nkx6-2" "Rgs1" "Pcdh10" "Zg16"
## [187] "Olfm4" "Ptcra" "Atp1b2" "Hmgn2" "Tulp3" "Gpr83"
## [193] "Klre1" "Cdkn1a" "Irf7" "Cebpd" "Timp3" "S1pr1"
## [199] "Lmna" "Peg3" "Ccr3" "Emp3" "Rnd3" "Aldh1a3"
## [205] "Foxf1" "H1f0" "Rgs16" "Gla" "Pfn1" "Kif4"
## [211] "Sox4" "Nmrk1" "Palld" "Coro1a" "Aspm" "H2afx"
## [217] "Kcnq1ot1" "Dusp2" "Plk2" "Plek" "Asb4" "Bst2"
## [223] "Gbp2" "Hoxa7" "Tshz2" "Ccr5" "Asf1b" "Clu"
## [229] "Foxd1" "Fzd8" "Tnfaip2" "Marcksl1" "Rep15" "Cd79a"
## [235] "Zfp503" "Lig1" "S100a10" "Foxp3" "Ccrl2" "Tk1"
## [241] "Serpinb2" "Irf4" "Klrg1" "Csn3" "Ms4a6b" "Rgs6"
## [247] "Cd40lg" "Ccr1" "Tsc22d1" "Pmepa1" "Bmp2" "Gimap7"
## [253] "Gbp2b" "Dapl1" "Tmem132c" "Pbk" "Npy" "Tchh"
## [259] "Pcdh18" "Elavl2" "Fam13a" "Ezh2" "Tmem158" "Runx1t1"
## [265] "Ar" "Cxcr6" "Cwc25" "Ptma" "Ccdc34" "Klra9"
## [271] "Eomes" "Gzmk" "Rgs10" "Vgf" "Serpina3f" "Adamts1"
## [277] "Irf8" "Rilpl2" "Anp32b" "Tagln" "Rab27b" "Tubb2b"
## [283] "Lat" "Anxa2" "Ndrg1" "Cd6" "Gpx1" "Nr4a1"
## [289] "Got1" "Lmnb1" "Cidea" "Epha5" "Tnf" "Klri2"
## [295] "Gbp6" "Cenpm" "Olfr889" "Hmgn3" "Klf4" "Sdc4"DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)
well, there’s nearly no big difference in PC1/2, a little bit shift in PC3/4, comparing with initial run
DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)
decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)
PCs <- 1:27
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 10)## Computing nearest neighbor graph## Computing SNNGEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 2.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 22088
## Number of edges: 306841
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7802
## Number of communities: 31
## Elapsed time: 1 secondsRun UMAP/tSNE
GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 998)## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 15:27:31 UMAP embedding parameters a = 0.9922 b = 1.112## 15:27:31 Read 22088 rows and found 27 numeric columns## 15:27:31 Using Annoy for neighbor search, n_neighbors = 20## 15:27:31 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 15:27:34 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\Rtmp2phpoV\file95fc4de13e2a
## 15:27:34 Searching Annoy index using 1 thread, search_k = 2000
## 15:27:39 Annoy recall = 100%
## 15:27:40 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 15:27:41 Initializing from normalized Laplacian + noise (using irlba)
## 15:27:41 Commencing optimization for 200 epochs, with 655162 positive edges
## 15:28:04 Optimization finishedDimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno", label.size = 2.75) +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno", label.size = 2.75) +
DimPlot(GEX.seur, reduction = "umap", label = T)
sort and test
GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
levels = c(1,8, # Tn
27, 19, # Treg
16, # Th17
25, # Th.CC
17, # Th1
30, # Th2
29, # ?
7, # iNKT
21, # MAIT
20,5, # gdT
26,10,28, # ILC2
24, # ILC3.CC
14,22,18,2,4,6,9,23, # ILC3.Ncr1
0, # ILC3.nn
3,13,12,15,11 # ILC3.Ccr6
))VlnPlot(subset(GEX.seur, subset = cnt %in% c("CTL")),
features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")
VlnPlot(subset(GEX.seur, subset = cnt %in% c("CKO")),
features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")
markers
# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "sort_clusters"
#GEX.markers.sort <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.05,
# test.use = "MAST",
# #test.use = "wilcox",
# logfc.threshold = 0.25)
GEX.markers.sort <- read.table("20230927_10x_SZJ.t1_adv.sort_markers.csv", header = TRUE, sep = ",")
#GEX.markers.sort %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)GEX.markers.sort$cluster <- factor(as.character(GEX.markers.sort$cluster),
levels = levels(GEX.seur$sort_clusters))
markers.sort_t60 <- (GEX.markers.sort %>% group_by(cluster) %>%
filter(pct.1>0.01 & gene %in% grep("Rps|Rpl|mt-",GEX.markers.sort$gene,invert = T,value = T)) %>%
top_n(n = 60, wt = avg_log2FC) %>%
filter(p_val_adj < 0.01) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.sort_t120 <- (GEX.markers.sort %>% group_by(cluster) %>%
filter(pct.1>0.01 & gene %in% grep("Rps|Rpl|mt-",GEX.markers.sort$gene,invert = T,value = T)) %>%
top_n(n = 120, wt = avg_log2FC) %>%
filter(p_val_adj < 0.01) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene747/60## [1] 12.45747/64## [1] 11.67188747/65## [1] 11.49231pp.t60 <- list()
for(i in 1:12){
pp.t60[[i]] <- DotPlot(GEX.seur, features = rev(markers.sort_t60[(64*i-63):(64*i)])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}
pp.t60## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
nnn = 1278
nnn/60## [1] 21.3nnn/64## [1] 19.96875nnn/65## [1] 19.66154pp.t120 <- list()
for(i in 1:20){
pp.t120[[i]] <- DotPlot(GEX.seur, features = rev(markers.sort_t120[(64*i-63):(64*i)])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}## Warning in FetchData.Seurat(object = object, vars = features, cells = cells):
## The following requested variables were not found: NA#pp.t120some reference papers for innate-like Tcell markers
Shilpi Chandra et al. Sci Immunol. 2023
https://pubmed.ncbi.nlm.nih.gov/37948512/
S Harsha Krovi NatCommun. 2020
https://pubmed.ncbi.nlm.nih.gov/33288744/
MAIT: Klrd1,Klra9,Nkg7,Ccl5,Cd7,Cd160,Xcl1 …
iNKT: Klrb1c,Ifng,Tbx21,Gzma …
still, here iNKT/MAIT clusters couldn’t be 100% confirmed,
because they usually are defined by specific TCR v-j genes
VlnPlot(GEX.seur, features = "Pdcd1", group.by = "sort_clusters") + NoLegend()
VlnPlot(GEX.seur, features = c("Fgl2","Klrd1","Klrk1","Klra9",
"Styk1","Nkg7","Ccl5","Cd7"), group.by = "sort_clusters",ncol = 1) + NoLegend()
VlnPlot(GEX.seur, features = c("Cd160","Klrb1c","Xcl1","Cxcr3","Gimap4"), group.by = "sort_clusters", ncol = 1) + NoLegend()
VlnPlot(GEX.seur, features = c("Klrb1c","Ifng","Tbx21","Gzma"), group.by = "sort_clusters",ncol = 1) + NoLegend()
VlnPlot(GEX.seur, features = "Eomes", group.by = "sort_clusters") + NoLegend()
DotPlot(GEX.seur, features = markers.new, group.by = "sort_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) +
scale_y_discrete(limits=rev)
Anno
# remove very few low-quality mix-like C29
GEX.seur <- subset(GEX.seur, subset = seurat_clusters != 29)
GEX.seur## An object of class Seurat
## 16621 features across 22038 samples within 1 assay
## Active assay: RNA (16621 features, 1197 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapGEX.seur@meta.data[,grep("pANN|snn",colnames(GEX.seur@meta.data),value = T)] <- NULL
head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAACAGGA-1 ILC3_GEX 3483 1341 2.871088 29.08412
## AAACCCAAGAATTGCA-1 ILC3_GEX 2244 1068 4.099822 19.51872
## AAACCCAAGCAGCCCT-1 ILC3_GEX 2074 981 3.905497 23.09547
## AAACCCAAGCCTGAAG-1 ILC3_GEX 3232 1370 2.165842 22.92698
## AAACCCAAGCGCATCC-1 ILC3_GEX 2902 1288 2.446589 23.56995
## AAACCCACAGAAGTGC-1 ILC3_GEX 2306 1142 2.428448 21.50911
## FB.info S.Score G2M.Score Phase seurat_clusters cnt
## AAACCCAAGAACAGGA-1 CKO.2 0.009665415 -0.04022873 S 19 CKO
## AAACCCAAGAATTGCA-1 CTL.3 -0.050866345 -0.08275860 G1 6 CTL
## AAACCCAAGCAGCCCT-1 CKO.4 -0.069848995 -0.07870466 G1 0 CKO
## AAACCCAAGCCTGAAG-1 CKO.1 -0.040119167 -0.10621201 G1 20 CKO
## AAACCCAAGCGCATCC-1 CTL.2 -0.024552021 -0.03206281 G1 20 CTL
## AAACCCACAGAAGTGC-1 CTL.3 -0.068676919 -0.11434633 G1 13 CTL
## DoubletFinder0.05 DoubletFinder0.1 sort_clusters preAnno
## AAACCCAAGAACAGGA-1 Singlet Singlet 19 Treg
## AAACCCAAGAATTGCA-1 Singlet Singlet 6 ILC3.Ncr1
## AAACCCAAGCAGCCCT-1 Singlet Singlet 0 ILC3.nn
## AAACCCAAGCCTGAAG-1 Singlet Singlet 20 gdT
## AAACCCAAGCGCATCC-1 Singlet Singlet 20 gdT
## AAACCCACAGAAGTGC-1 Singlet Singlet 13 ILC3.Ccr6GEX.seur$Anno <- as.character(GEX.seur$seurat_clusters)
## Cd3d+Cd4+
GEX.seur$Anno[GEX.seur$Anno %in% c(1,8)] <- "Tn" # Cd4+(weak)Sell+Ccr7+
GEX.seur$Anno[GEX.seur$Anno %in% c(27,19)] <- "Treg" # Foxp3+
GEX.seur$Anno[GEX.seur$Anno %in% c(16)] <- "Th17" # Il17a+
GEX.seur$Anno[GEX.seur$Anno %in% c(25)] <- "Th.CC" # cycling, most should be Th17
GEX.seur$Anno[GEX.seur$Anno %in% c(17)] <- "Th1" # Gzmk+
GEX.seur$Anno[GEX.seur$Anno %in% c(30)] <- "Th2" # Cd3+Cd4+ and co-express many ILC2 markers
## Cd3d+Cd4-
GEX.seur$Anno[GEX.seur$Anno %in% c(7)] <- "iNKT" # Cd160+Tbx21+Klrb1c+
GEX.seur$Anno[GEX.seur$Anno %in% c(21)] <- "MAIT" # Cd160+Klrd1+Cd7+Itgae+
GEX.seur$Anno[GEX.seur$Anno %in% c(20,5)] <- "gdT" # Tcrg-C1+Cd163l1+Trdv4+, C20 Il17a+
## Cd3d-
GEX.seur$Anno[GEX.seur$Anno %in% c(26,10,28)] <- "ILC2" # Gata3+Calca+Il5+
# Rorc+
GEX.seur$Anno[GEX.seur$Anno %in% c(24)] <- "ILC3.CC" # cycling, should be Ncr1+
GEX.seur$Anno[GEX.seur$Anno %in% c(14,22,18,2,4,6,9,23)] <- "ILC3.Ncr1" # Ncr1+
GEX.seur$Anno[GEX.seur$Anno %in% c(0,3)] <- "ILC3.nn" # Ncr1-Ccr6- double negative
GEX.seur$Anno[GEX.seur$Anno %in% c(13,12,15,11)] <- "ILC3.Ccr6" # Ccr6+
GEX.seur$Anno <- factor(GEX.seur$Anno,
levels = c("Tn","Treg","Th17","Th.CC","Th1","Th2",
"iNKT","MAIT","gdT",
"ILC2",
"ILC3.CC","ILC3.Ncr1","ILC3.nn","ILC3.Ccr6"))DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno", label.size = 2.75) +
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters")
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno", label.size = 2.75) +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
Idents(GEX.seur) <- "Anno"
GEX.seur$rep <- paste0("rep",
gsub("CTL|CKO","",as.character(GEX.seur$FB.info)))
#saveRDS(GEX.seur, "./Baf53cre-AhrCKO.SI_ILC3CD4T.Anno.rds")final markers
## define feature colors
# Cell2020
material.heat = function(n){
colorRampPalette(
c(
"#283593", # indigo 800
"#3F51B5", # indigo
"#2196F3", # blue
"#00BCD4", # cyan
"#4CAF50", # green
"#8BC34A", # light green
"#CDDC39", # lime
"#FFEB3B", # yellow
"#FFC107", # amber
"#FF9800", # orange
"#FF5722" # deep orange
)
)(n)
}pp.Anno <- DotPlot(GEX.seur, features = markers.new, group.by = "Anno",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) +
scale_color_gradientn(colours = material.heat(100)) +
scale_y_discrete(limits=rev)
pp.Anno
markers.final <- c("Cd3d","Cd4","Ccr7","Sell",
"Hopx","Foxp3","Il10",
"Il17a","Nebl","Top2a","Mki67",
"Gzma","Gzmk","Ccr5","Il12rb2","Ifng",
"Tbx21","Klrb1c","Xcl1","Nkg7","Plac8","Ccl5",
"Cd160","Klrd1","Cd7","Itgae",
"Tcrg-C1","Cd163l1","Blk",
"Gata3","Calca","Il5","Il13",
"Ncr1","Gzmb","Klrk1","Fcer1g",
"Il22","Ccr6","Rorc","Il17f")pn.Anno.a <- DotPlot(GEX.seur, features = markers.final, group.by = "Anno",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) +
scale_color_gradientn(colours = material.heat(100)) +
scale_y_discrete(limits=rev)
pn.Anno.a
pn.Anno.b <- DotPlot(GEX.seur, features = markers.final, group.by = "Anno",
cols = c("midnightblue","darkorange1")) +
coord_flip() +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1, size = 9)) +
scale_color_gradientn(colours = material.heat(100)) #+
#scale_y_discrete(limits=rev)
pn.Anno.b
DEGs
#
Idents(GEX.seur) <- "cnt"all.clusters <- levels(GEX.seur$Anno)
all.clusters## [1] "Tn" "Treg" "Th17" "Th.CC" "Th1" "Th2"
## [7] "iNKT" "MAIT" "gdT" "ILC2" "ILC3.CC" "ILC3.Ncr1"
## [13] "ILC3.nn" "ILC3.Ccr6"#
list_new <- list()
list_new[["All"]] <- all.clusters
list_new[["CD4T"]] <- all.clusters[1:6]
list_new[["ILC3"]] <- all.clusters[11:14]
for(NN in all.clusters){
list_new[[NN]] <- NN
}
names_new <- names(list_new)
list_new## $All
## [1] "Tn" "Treg" "Th17" "Th.CC" "Th1" "Th2"
## [7] "iNKT" "MAIT" "gdT" "ILC2" "ILC3.CC" "ILC3.Ncr1"
## [13] "ILC3.nn" "ILC3.Ccr6"
##
## $CD4T
## [1] "Tn" "Treg" "Th17" "Th.CC" "Th1" "Th2"
##
## $ILC3
## [1] "ILC3.CC" "ILC3.Ncr1" "ILC3.nn" "ILC3.Ccr6"
##
## $Tn
## [1] "Tn"
##
## $Treg
## [1] "Treg"
##
## $Th17
## [1] "Th17"
##
## $Th.CC
## [1] "Th.CC"
##
## $Th1
## [1] "Th1"
##
## $Th2
## [1] "Th2"
##
## $iNKT
## [1] "iNKT"
##
## $MAIT
## [1] "MAIT"
##
## $gdT
## [1] "gdT"
##
## $ILC2
## [1] "ILC2"
##
## $ILC3.CC
## [1] "ILC3.CC"
##
## $ILC3.Ncr1
## [1] "ILC3.Ncr1"
##
## $ILC3.nn
## [1] "ILC3.nn"
##
## $ILC3.Ccr6
## [1] "ILC3.Ccr6"names_new## [1] "All" "CD4T" "ILC3" "Tn" "Treg" "Th17"
## [7] "Th.CC" "Th1" "Th2" "iNKT" "MAIT" "gdT"
## [13] "ILC2" "ILC3.CC" "ILC3.Ncr1" "ILC3.nn" "ILC3.Ccr6"run
#test.DEGs_new# save DEGs' table
df_test.DEGs_new <- data.frame()
for(nn in names_new){
df_test.DEGs_new <- rbind(df_test.DEGs_new, data.frame(test.DEGs_new[[nn]],Anno=nn))
}
#write.csv(df_test.DEGs_new, "Baf53cre_AhrCKO.SI_ILC3CD4T.DEGs_CTLvsCKO_Anno.csv")stat
cut0: padj 0.05 log2FC 0.1
cut1: padj 0.05 log2FC 0.2
cut2: padj 0.01 log2FC > log2(1.5) (|FC|>1.5)
cut0
# cut0
cut.padj = 0.05
cut.log2FC = 0.1
cut.pct1 = 0.1stat_test0.DEGs_new <- df_test.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(Anno,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
stat_test0.DEGs_new## Anno cluster up.DEGs
## 1 All CTL 50
## 2 All CKO 19
## 3 CD4T CTL 4
## 4 CD4T CKO 8
## 5 gdT CKO 4
## 6 ILC2 CKO 1
## 7 ILC3 CTL 39
## 8 ILC3 CKO 19
## 9 ILC3.Ccr6 CKO 3
## 10 ILC3.Ncr1 CTL 26
## 11 ILC3.Ncr1 CKO 19
## 12 ILC3.nn CTL 3
## 13 ILC3.nn CKO 3
## 14 iNKT CKO 1
## 15 Th1 CTL 1
## 16 Tn CTL 1
## 17 Tn CKO 3df_test.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(Anno,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=Anno, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt, name="") +
labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
cut1
# cut1
cut.padj = 0.05
cut.log2FC = 0.2
cut.pct1 = 0.1stat_test1.DEGs_new <- df_test.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(Anno,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
stat_test1.DEGs_new## Anno cluster up.DEGs
## 1 All CTL 4
## 2 All CKO 3
## 3 CD4T CTL 3
## 4 CD4T CKO 6
## 5 gdT CKO 2
## 6 ILC2 CKO 1
## 7 ILC3 CTL 6
## 8 ILC3 CKO 3
## 9 ILC3.Ccr6 CKO 2
## 10 ILC3.Ncr1 CTL 5
## 11 ILC3.Ncr1 CKO 5
## 12 ILC3.nn CTL 3
## 13 ILC3.nn CKO 1
## 14 iNKT CKO 1
## 15 Th1 CTL 1
## 16 Tn CTL 1
## 17 Tn CKO 2df_test.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(Anno,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=Anno, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt, name="") +
labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
plot
pp_test.DEGs.new <- lapply(list_new, function(x){NA})#
test.DEGs_new.combine <- test.DEGs_new
test.DEGs_new.combine <- lapply(test.DEGs_new.combine, function(x){
x[x$cluster=="CTL","avg_log2FC"] <- -x[x$cluster=="CTL","avg_log2FC"]
x
})All
pp_test.DEGs.new$All <- test.DEGs_new.combine$All %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="All CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$All
CD4T
pp_test.DEGs.new$CD4T <- test.DEGs_new.combine$CD4T %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="CD4T CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$CD4T
ILC3
pp_test.DEGs.new$ILC3 <- test.DEGs_new.combine$ILC3 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="ILC3 CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3
Tn
pp_test.DEGs.new$Tn <- test.DEGs_new.combine$Tn %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="Tn CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Tn
Treg
pp_test.DEGs.new$Treg <- test.DEGs_new.combine$Treg %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="Treg CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Treg
Th17
pp_test.DEGs.new$Th17 <- test.DEGs_new.combine$Th17 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="Th17 CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Th17
Th.CC
pp_test.DEGs.new$Th.CC <- test.DEGs_new.combine$Th.CC %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="Th.CC CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Th.CC
Th1
pp_test.DEGs.new$Th1 <- test.DEGs_new.combine$Th1 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="Th1 CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Th1
Th2
pp_test.DEGs.new$Th2 <- test.DEGs_new.combine$Th2 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="Th2 CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$Th2
iNKT
pp_test.DEGs.new$iNKT <- test.DEGs_new.combine$iNKT %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="iNKT CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$iNKT
MAIT
pp_test.DEGs.new$MAIT <- test.DEGs_new.combine$MAIT %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="MAIT CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$MAIT
gdT
pp_test.DEGs.new$gdT <- test.DEGs_new.combine$gdT %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="gdT CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$gdT
ILC2
pp_test.DEGs.new$ILC2 <- test.DEGs_new.combine$ILC2 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="ILC2 CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC2
ILC2
pp_test.DEGs.new$ILC2 <- test.DEGs_new.combine$ILC2 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="ILC2 CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC2
ILC3.CC
pp_test.DEGs.new$ILC3.CC <- test.DEGs_new.combine$ILC3.CC %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="ILC3.CC CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3.CC
ILC3.Ncr1
pp_test.DEGs.new$ILC3.Ncr1 <- test.DEGs_new.combine$ILC3.Ncr1 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="ILC3.Ncr1 CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3.Ncr1
ILC3.nn
pp_test.DEGs.new$ILC3.nn <- test.DEGs_new.combine$ILC3.nn %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="ILC3.nn CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3.nn
ILC3.Ccr6
pp_test.DEGs.new$ILC3.Ccr6 <- test.DEGs_new.combine$ILC3.Ccr6 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.1)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Gm",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("CTL"=as.vector(color.cnt[1]),
"CKO"=as.vector(color.cnt[2]),
"None"="grey")) +
theme_classic() + labs(title="ILC3.Ccr6 CTL vs CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test.DEGs.new$ILC3.Ccr6
composition
all
stat1_Anno <- GEX.seur@meta.data[,c("cnt","rep","Anno")]
stat1_Anno.s <- stat1_Anno %>%
group_by(cnt,rep,Anno) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.pcom1.Anno <- stat1_Anno.s %>%
ggplot(aes(x = Anno, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.78, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt, name="") +
labs(title="Cell Composition of SI_ILC3CD4T.CTL_CKO, Anno", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=Anno, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom1.Anno
sig.test
glm.nb - anova.LRT
Sort.N <- c("CTL","CKO")
glm.nb_anova.1.Anno <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat1_Anno.s_N <- stat1_Anno.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat1_Anno.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat1_Anno.s_N$total <- total.N[paste0(stat1_Anno.s_N$cnt,
"_",
stat1_Anno.s_N$rep),"total"]
glm.nb_anova.1.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat1_Anno.s_N$Anno)){
glm.nb_anova.1.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat1_Anno.s_N %>% filter(Anno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat1_Anno.s_N %>% filter(Anno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat1_Anno.s_N %>% filter(Anno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.1.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.1.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}glm.nb_anova.1.Anno_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.1.Anno)))
rownames(glm.nb_anova.1.Anno_df) <- names(glm.nb_anova.1.Anno)
colnames(glm.nb_anova.1.Anno_df) <- gsub("X","C",colnames(glm.nb_anova.1.Anno_df))
glm.nb_anova.1.Anno_df## Tn Treg Th17 Th.CC Th1 Th2
## CTLvsCKO 0.3912709 0.008942821 0.8056086 0.8006229 0.7828195 0.03982599
## iNKT MAIT gdT ILC2 ILC3.CC ILC3.Ncr1
## CTLvsCKO 0.7431366 0.6684469 0.01670645 0.5941181 0.01892623 0.6151186
## ILC3.nn ILC3.Ccr6
## CTLvsCKO 0.8142205 0.7319146round(glm.nb_anova.1.Anno_df,4)## Tn Treg Th17 Th.CC Th1 Th2 iNKT MAIT gdT ILC2
## CTLvsCKO 0.3913 0.0089 0.8056 0.8006 0.7828 0.0398 0.7431 0.6684 0.0167 0.5941
## ILC3.CC ILC3.Ncr1 ILC3.nn ILC3.Ccr6
## CTLvsCKO 0.0189 0.6151 0.8142 0.7319CD4T
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno", label.size = 2.75,
cols = c(scales::hue_pal()(14)[c(1:6)],rep("lightgrey",8) )) +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
list_new$CD4T## [1] "Tn" "Treg" "Th17" "Th.CC" "Th1" "Th2"stat2_Anno <- GEX.seur@meta.data[,c("cnt","rep","Anno")] %>%
filter(Anno %in% list_new$CD4T)
stat2_Anno.s <- stat2_Anno %>%
group_by(cnt,rep,Anno) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.pcom2.Anno <- stat2_Anno.s %>%
ggplot(aes(x = Anno, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.78, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt, name="") +
labs(title="Cell Composition, CD4T.CTL_CKO, Anno", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=Anno, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom2.Anno
sig.test
glm.nb - anova.LRT
Sort.N <- c("CTL","CKO")
glm.nb_anova.2.Anno <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat2_Anno.s_N <- stat2_Anno.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat2_Anno.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat2_Anno.s_N$total <- total.N[paste0(stat2_Anno.s_N$cnt,
"_",
stat2_Anno.s_N$rep),"total"]
glm.nb_anova.2.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat2_Anno.s_N$Anno)){
glm.nb_anova.2.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat2_Anno.s_N %>% filter(Anno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat2_Anno.s_N %>% filter(Anno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat2_Anno.s_N %>% filter(Anno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.2.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.2.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}glm.nb_anova.2.Anno_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.2.Anno)))
rownames(glm.nb_anova.2.Anno_df) <- names(glm.nb_anova.2.Anno)
colnames(glm.nb_anova.2.Anno_df) <- gsub("X","C",colnames(glm.nb_anova.2.Anno_df))
glm.nb_anova.2.Anno_df## Tn Treg Th17 Th.CC Th1 Th2 iNKT
## CTLvsCKO 0.736032 0.0006323642 0.9614593 0.3357168 0.1291926 0.002086836 NA
## MAIT gdT ILC2 ILC3.CC ILC3.Ncr1 ILC3.nn ILC3.Ccr6
## CTLvsCKO NA NA NA NA NA NA NAround(glm.nb_anova.2.Anno_df,4)## Tn Treg Th17 Th.CC Th1 Th2 iNKT MAIT gdT ILC2 ILC3.CC
## CTLvsCKO 0.736 6e-04 0.9615 0.3357 0.1292 0.0021 NA NA NA NA NA
## ILC3.Ncr1 ILC3.nn ILC3.Ccr6
## CTLvsCKO NA NA NAILCs
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "Anno", label.size = 2.75,
cols = c(rep("lightgrey",9) ,scales::hue_pal()(14)[c(10:14)])) +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
list_new$ILC3## [1] "ILC3.CC" "ILC3.Ncr1" "ILC3.nn" "ILC3.Ccr6"stat3_Anno <- GEX.seur@meta.data[,c("cnt","rep","Anno")] %>%
filter(Anno %in% c(list_new$ILC3,"ILC2"))
stat3_Anno.s <- stat3_Anno %>%
group_by(cnt,rep,Anno) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.pcom3.Anno <- stat3_Anno.s %>%
ggplot(aes(x = Anno, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.78, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt, name="") +
labs(title="Cell Composition, ILCs.CTL_CKO, Anno", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=Anno, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom3.Anno
sig.test
glm.nb - anova.LRT
Sort.N <- c("CTL","CKO")
glm.nb_anova.3.Anno <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat3_Anno.s_N <- stat3_Anno.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat3_Anno.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat3_Anno.s_N$total <- total.N[paste0(stat3_Anno.s_N$cnt,
"_",
stat3_Anno.s_N$rep),"total"]
glm.nb_anova.3.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat3_Anno.s_N$Anno)){
glm.nb_anova.3.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat3_Anno.s_N %>% filter(Anno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat3_Anno.s_N %>% filter(Anno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat3_Anno.s_N %>% filter(Anno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.3.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.3.Anno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}glm.nb_anova.3.Anno_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.3.Anno)))
rownames(glm.nb_anova.3.Anno_df) <- names(glm.nb_anova.3.Anno)
colnames(glm.nb_anova.3.Anno_df) <- gsub("X","C",colnames(glm.nb_anova.3.Anno_df))
glm.nb_anova.3.Anno_df## Tn Treg Th17 Th.CC Th1 Th2 iNKT MAIT gdT ILC2 ILC3.CC
## CTLvsCKO NA NA NA NA NA NA NA NA NA 0.4784894 0.01651209
## ILC3.Ncr1 ILC3.nn ILC3.Ccr6
## CTLvsCKO 0.3515031 0.6700376 0.7037525round(glm.nb_anova.3.Anno_df,4)[,10:14,drop=F]## ILC2 ILC3.CC ILC3.Ncr1 ILC3.nn ILC3.Ccr6
## CTLvsCKO 0.4785 0.0165 0.3515 0.67 0.7038